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We present the application of a new method to compute the Wiener filter solution of large and 
complex data sets. Contrary to the iterative solvers usually employed in signal processing, our 



iy) algorithm does not require the use of preconditioners to be computationally efficient. The new 

i scheme is conceptually very simple and therefore easy to implement, numerically absolutely sta- 

_ . ble, and guaranteed to converge. We introduce a messenger field to mediate between the different 

^ preferred bases in which signal and noise properties can be specified most conveniently, and 

rephrase the signal reconstruction problem in terms of this auxiliary variable. We demonstrate the 
capabilities of the algorithm by applying it to cosmic microwave background (CMB) radiation 
data obtained by the WMAP satellite. 
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1. Introduction 

Among the linear filters that make use of statistical information about the data, the generalized 
Wiener filter stands out as the maximum a posteriori solution for the case that signal and noise 
both follow a Gaussian distribution (Wiener 1949). Although Wiener filtering is of fundamental 
importance, only recently efficient algorithms have been developed for the relaxed assumption that 
the signal covariance is not known in advance (e.g., Wandelt et al. 2004; EnBlin & Frommert 201 1). 

Assuming a linear model for the observed data d as a combination of signal s and noise n, 

d — s + n^ (1.1) 

the Wiener filter swf is defined as the solution of the equation 

(S^+N-^wf^N- 1 ;/. (1.2) 

It plays an important role in signal processing. Taking the example of cosmic microwave back- 
ground (CMB) radiation data, it is calculated during optimal power spectrum estimation (e.g., 
Tegmark 1997b; Bond et al. 1998; Oh et al. 1999; Eisner & Wandelt 2012b), likelihood analy- 
sis (e.g., Hinshaw et al. 2007; Dunkley et al. 2009; Eisner & Wandelt 2012c), mapmaking (e.g., 
Tegmark et al. 1997; Tegmark 1997a), and lensing reconstructions (e.g., Hirata et al. 2004; Smith 
et al. 2007), etc. 

To evaluate Eq. (1.2) in practice can turn out very challenging for the large and complex data 
sets obtained by state of the art experiments. Problems occur as the size of S and N increase with 
the second power of the number of data samples, rendering the storage and processing of dense 
systems impractical. Fortunately, it is often feasible to specify sets of bases where signal and 
noise covariance become sparse. However, these bases are generally incompatible, i.e., it is not 
possible to represent them in a single basis as sparse systems such that Eq. (1.2) can be solved 
trivially. Numerical algorithms for the exact solution of the Wiener filter equation are therefore 
often complex, e.g., involving conjugate gradient solvers with multigrid preconditioners (Smith 
et al. 2007). Finding fast preconditioners that work is a highly non-trivial art, especially since the 
matrices of interest are often extremely ill-conditioned (Kn^/Knin ^ 10 7 ). 

In Sect. 2, we briefly describe a new iterative algorithm for the solution of the exact Wiener 
filter equation. We illustrate the approach by applying it to WMAP7 data in Sect. 3. Summarizing 
the main aspects of our findings, we then discuss possible extensions and areas of applications 
(Sect. 4). Details on the mathematical aspects of the method can be found in Eisner & Wandelt 
(2013). 

2. Method 

To efficiently compute the Wiener filter, we first introduce a messenger field t with covariance 
T. The covariance properties of this auxiliary variable are very simple: T is proportional to the 
identity matrix, a property that is preserved under any orthogonal basis change. Therefore, while 
it may not be possible to directly apply expressions like (S + N) _1 to a data vector, we can always 
do so with combinations as (S + T) _1 and (N + T) _1 , no matter what basis is chosen to render S 
and N sparse. 
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To benefit from this possibility, we now specify equations for the signal reconstruction s and 
the auxiliary field t, 



where we defined N = N — T and introduced one additional scalar parameter X, which we will 
use to accelerate convergence. We specify the covariance matrix of the auxiliary field according 
to T = min(diag(N)) • 1. In the limit X = 1, the system of equations defined above reduces to the 
Wiener filter equation. 

An outline of the algorithm can be given as follows. Initially, the vectors s and t are set to 
zero. We solve Eq. (2.1) for the auxiliary field t in the basis defined by the noise covariance matrix. 
Next, we change the basis to, e.g., Fourier space, where S can be represented easily. Then, we 
solve for s using Eq. (2.2) given the latest realization of t, and finally transform the result back 
to the original basis. It can be shown that the signal reconstruction s converges exponentially and 
unconditionally to the Wiener filter solution swf (Eisner & Wandelt 2013). Based on a comparison 
to more standard conjugate gradient solvers, we find the final map to be accurate to about 1 part in 
10 5 , depending on the adopted stopping criterion. 

3. WMAP analysis 

We now demonstrate the application of the algorithm to the V-band data of the Wilkinson mi- 
crowave anisotropy probe (WMAP7 1 , Jarosik et al. 201 1). We adopted the official WMAP extended 
temperature and polarization masks, reducing the sky fraction to around 70 %. 

For the best performance in terms of computing time, we changed X in the following way: 
To initialize the algorithm, we chose it according to the ratio T/Cj^ at £ litx = 20. As shown in 
Fig. 1, a high value of X corresponds to a large convolution kernel in Eq. (2.2). Accordingly, small 
scales are smoothed out and information is more efficiently propagated into the masked regions. To 
monitor convergence of the algorithm, we used the % 2 = s^S~ l s + (d — s)^N _1 (d — s) of the current 
solution s. After it has reached a plateau, we decreased X to the ratio at £f™ = 2£f^ r The evolution 
of X, together with the current value of % 2 at every iteration, is depicted in Fig. 2. We stopped 
the algorithm as soon as the quality of the signal reconstruction improved only marginally between 
iterations, A% 2 < 10~ 4 0^2 at X = 1, and obtain a final goodness of fit of X 2 / n d.o.t = 0.996. 

Since we can attribute a specific length scale to a given value of X (the current value of AterX 
and fluctuations on smaller scales are suppressed owing to the width of the convolution kernel, we 
can restrict the maximum multipole moment of the computationally expensive spherical harmonic 
transforms to £^J « £ litx + 100. As a result, the overall computing time for the spherical harmonic 
transforms can be reduced by a large margin. We note that for the most expensive final iterations, 
where £^J « £ m2iX and more than half of the computing time is spent, the algorithm is an excellent 
candidate for further acceleration by means of graphics processing units (Eisner & Wandelt 2011). 

The buildup of the WMAP temperature Wiener filtered map is shown in Fig. 3. Characteristic 
for Wiener filter solutions, large scale fluctuations extend well inside masked regions. In Fig. 4, 

Available from http : / /lambda . gsf c . nasa . gov 



(N _1 + (AT)" 1 ) t = N" 1 d+ (AT)" 1 s 
(S- 1 + (lT)- 1 )s=(XT)- 1 t, 



(2.1) 
(2.2) 
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Figure 1: Convolution kernel. For A = 1 {blue lines) and A = 1000 {black lines), we show the shape of 
the convolution kernel used to solve for the signal given the current reconstruction of the messenger field in 
harmonic space {left panel) and real space {right panel). 




Iteration Iteration 



Figure 2: Convergence of the algorithm. Left panel: Starting from ^(lO 4 ), we reduce A step by step to 
unity. Right panel: For a given A , the X 2 of the solution first drops quickly and then reaches a plateau {blue 
line). The expectation value of X 2 for the final solution, given by the number of degrees of freedom, is also 
indicated {black line). 

we plot the power spectra of the reconstruction for temperature and polarization, and demonstrate 
that the Wiener filter maps can be augmented to constrained realizations, with the correct signal 
variances. 

4. Conclusion 

We have summarized a new method to efficiently calculate the Wiener filter solution of general 
data sets. As a sample application, we analyzed WMAP temperature and polarization maps. 

The algorithm is not only simple to implement, but also robust. Even after including polar- 
ization data in the analysis, it maintains good performance, despite of the increase in the condition 
number of the covariance matrices. We consider this fact to be of particular importance as conju- 
gate gradient solvers are likely to perform noticeably worse in this setup. 
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Figure 3: Buildup of the temperature reconstruction. Snapshots of the algorithm after 50 iterations {upper 
panel), 250 iterations {middle panel), and of the final solution {lower panel) show the convergence of the 
reconstructed signal from the largest to the smallest scales. 
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Figure 4: Power spectra of the reconstruction. In this series of plots of the temperature (TT, upper panel), 
polarization (EE, middle panel), and cross power spectra (TE, lower panel), we show that a constrained 
realization obtained with this algorithm (red line), consisting of the the sum of the Wiener filter solution (blue 
line) and a fluctuation map (green line), is unbiased compared to the fiducial power spectrum multiplied by 
the beam function (black line). 
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We finally note that the algorithm is also very flexible, and many possible extensions are still 
to be explored. For some problems, for example, it may be beneficial to solve for the noise vector 
instead of the signal, which is calculated indirectly from the result. The auxiliary field then becomes 
associated with the signal reconstruction instead of the noise. In certain situations, it may also prove 
useful to include more than one messenger field, for example if multiple observations (e.g., from 
different detectors) should be combined in a joint analysis. It also remains to be explored to what 
extent our method could be combined with conventional iterative schemes, for example to act as a 
smoother in a multigrid scheme to speed up convergence further. 
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